clear all
cls
global path "C:\...."
* dtapath contains the path to the folder with intermediate Stata data files
global dtapath "$path\dta-files"
*resulspath contain the path to the folder where all output from this do-file will be saved
global resultspath "$path\results"
capture mkdir "$resultspath"
cd "$resultspath"
version 16

global zeros="no" // use dataset with imputed zeros?


***********************************************************
*** SAME AGGREGATION LEVEL AS IN B&B BUT USING COMTRADE ***
***********  Replicates results in footnote 12  ***********
***********************************************************
cap log close

if "$zeros"=="yes" {
cap log close
log using BB_replication_withzeros, replace
use "$dtapath\BB_allyears_all_zeros_nomirror.dta", replace
}
else {
log using BB_replication_nozeros, replace
use "$dtapath\BB_allyears_all_nozeros_nomirror.dta", replace
}


	drop if year>2000  // B&B use data from 1960-2000. After we drop all years 
	* after 2000, we have data for 1962, 1965, 1970, 1975, 1980, 1985, 1990, and 2000.
	
	* collapse the data and take logs to get (log) aggregate bilateral trade flows (B&B's main dependent variable)
	collapse (sum) imports (mean)  comlang_off comcol curcol col45 distw EIA WTO_j WTO_i WTO PTA scott_fta fta, by(importeriso3 exporteriso3 year)
	
	egen pairid_str=concat(exporteriso3 importeriso3)
	encode pairid_str, gen(pairid)
	egen exp = concat(exporteriso3 year)
	egen imp = concat(importeriso3 year)
	encode exp, gen(ct_exp)
	encode imp, gen(ct_imp)
	xtset pairid year
	
	bys importeriso3 exporteriso3 (year): gen l1fta=fta[_n-1] if year==year[_n-1]+5 | (year==1965 & year==year[_n-1]+3)
	bys importeriso3 exporteriso3 (year): gen l2fta=fta[_n-2] if year==year[_n-2]+10 | (year==1970 & year==year[_n-2]+8)
	bys importeriso3 exporteriso3 (year): gen f1fta=fta[_n+1] if year==year[_n+1]-5 | (year==1962 & year==year[_n+1]-3)
	
	
	gen limports=log(imports)
	sum
	
	*  replicate key B&B regression (their Table 5, col. 4)
	
	reghdfe limports fta l1fta l2fta f1fta , abs(ct_exp ct_imp pairid)
	
	* total FTA effect 
	
	disp _b[fta]+_b[l1fta]+_b[l2fta]

	*for comparison, (B&B (2007, JIE) get a total FTA effect of 0.28+0.27+0.21 = 0.76 (log points)
	
log close
	

**************************************************************
*** DIFFERENT AGGREGATION LEVELS, HOMOGENEOUS COEFFICIENTS ***
*************** Replicates results in Table 1 ****************
**************************************************************
global zeros="yes" // use dataset with imputed zeros?

cd "$resultspath

capture log close
cls
disp("$zeros")

if "$zeros"=="yes" {

use "$dtapath\BB_allyears_all_zeros_nomirror.dta", replace
log using BB_FEcombis_withzeros, replace
}
else {

use "$dtapath\BB_allyears_all_nozeros_nomirror.dta", replace
log using BB_FEcombis_nozeros, replace
}


* drop all 4-digit sectors with less than 2,000 observations with positive trade flows;
* when there are too few positive flows, it is generally not possible to identify all parameters in the product-level regressions below 

drop if year>2000  // o/w counts not compatible with sector-by-sector regressions
gen sitc_4_digit=substr(commoditycode,4,4)
gen temp=(imports>0)
egen countpos=sum(temp),by(sitc_4_digit)
keep if countpos>=2000

/*
* additional check: drop  sectors that get dropped at various levels of aggregation 
* in the regressions below; goal is to keep the observations used at various levels
* comparable
* note: need to run the bit below first to generate the sample files
* aggregate data (SITC 0-digit: merge by country-pair and year since all 
* products are stored as "0" in the data at the zero digit level
cap drop _merge
merge n:1 exporteriso3 importeriso3 year using sample_0_digit.dta
keep if sample_0_digit==1
*/
* save under new name (use this file from now on!)
save "$dtapath\BB_allyears_all_zeros_2000obsplus.dta", replace

cd "$resultspath"
cls

global reglist="fta l1fta l2fta"

* 4-digit SITC: 625 products
foreach method in ols ppml {
*foreach method in ols {
	foreach fe in  imp_exp_bilateral all  {
		foreach level in 0 2 4 {	
	
	disp("`method'")
	disp("`fe'")
	disp("`level'")
	
			* collapse to right level
			use "$dtapath\BB_allyears_all_zeros_2000obsplus.dta", replace
			cap drop sitc*
			gen sitc_`level'_digit=substr(commoditycode,4,`level')
			replace sitc_`level'_digit="0" if sitc_`level'_digit==""
			drop if year>2000
			if `level'<4 {
			collapse (sum) imports (mean) fta, by(importeriso3 exporteriso3 sitc_`level'_digit year)
			}
						
			* generate lags and log imports
			if "`method'"=="ols" {
					*keep if imports~=0   // o/w we get more lags than in original B&B (2007) paper!
			}
			
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen l1fta=fta[_n-1] if year==year[_n-1]+5 | (year==1965 & year==year[_n-1]+3)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen l2fta=fta[_n-2] if year==year[_n-2]+10 | (year==1970 & year==year[_n-2]+8)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen f1fta=fta[_n+1] if year==year[_n+1]-5 | (year==1962 & year==year[_n+1]-3)
			gen limports=log(imports)
			sum limports imports fta l1fta l2fta 
			sum limports imports fta l1fta l2fta if limports~=.
			
			* ppml/ols regressions
			egen pairid_str=concat(exporteriso3 importeriso3)
			codebook pairid_str // check if number of country-pairs varies with aggregation
			encode pairid_str, gen(pairid)
			if "`fe'"=="none" {
				if "`method'"=="ppml" {
					poisson imports $reglist, cluster(pairid)
					
				}
				else if "`method'"=="ols" {
					reg limports $reglist, cluster(pairid)
				}	
				}
			
			if "`fe'"=="imp_exp" {
				egen exp_year=group(exporteriso3 year)
				egen imp_year=group(importeriso3 year)
				codebook exp_year imp_year
				if "`method'"=="ppml" {
				sdf
					ppmlhdfe imports $reglist, abs(exp_year imp_year) cluster(pairid) separation(fe mu)
				}
				else if "`method'"=="ols" {
					reghdfe limports $reglist, abs(exp_year imp_year) cluster(pairid)
				}
			}
			
			if "`fe'"=="imp_exp_bilateral" {
				egen exp_year=group(exporteriso3 year)
				egen imp_year=group(importeriso3 year)
				egen imp_exp=group(exporteriso3 importeriso3)
				codebook exp_year imp_year imp_exp
				if "`method'"=="ppml" {
					ppmlhdfe imports $reglist, abs(exp_year imp_year imp_exp) cluster(pairid) separation(fe mu)
				
				}
				else if "`method'"=="ols" {
					reghdfe limports $reglist, abs(exp_year imp_year imp_exp) cluster(pairid)
				}
			}
			
			
			if "`fe'"=="imp_exp_prod" {
				egen exp_year=group(exporteriso3 year)
				egen imp_year=group(importeriso3 year)
				egen prod_year=group(sitc_`level'_digit year)
				codebook exp_year imp_year prod_year
				if "`method'"=="ppml" {
					ppmlhdfe imports $reglist, abs(exp_year imp_year prod_year) cluster(pairid) separation(fe mu)
				}
				else if "`method'"=="ols" {
					reghdfe limports $reglist, abs(exp_year imp_year prod_year) cluster(pairid)
				}
			}
			if "`fe'"=="iprod_eprod" {
				egen exp_prod_year=group(exporteriso3 sitc_`level'_digit year)
				egen imp_prod_year=group(importeriso3 sitc_`level'_digit year)
				codebook exp_prod_year imp_prod_year
				if "`method'"=="ppml" {
					ppmlhdfe imports $reglist, abs(exp_prod_year imp_prod_year) cluster(pairid) separation(fe mu)
				}
				else if "`method'"=="ols" {
					reghdfe limports $reglist, abs(exp_prod_year imp_prod_year) cluster(pairid)
				}
			}
			if "`fe'"=="all" {
				egen exp_prod_year=group(exporteriso3 sitc_`level'_digit year)
				egen imp_prod_year=group(importeriso3 sitc_`level'_digit year)
				egen imp_exp_prod=group(importeriso3 exporteriso3 sitc_`level'_digit)
				codebook exp_prod_year imp_prod_year imp_exp_prod
				if "`method'"=="ppml" {
					ppmlhdfe imports $reglist, abs(exp_prod_year imp_prod_year imp_exp_prod) cluster(pairid) separation(fe mu)
				}
				else if "`method'"=="ols" {
					reghdfe limports $reglist, abs(exp_prod_year imp_prod_year imp_exp_prod) cluster(pairid)
				}
			}

			
				* total FTA effect and SE
				global tot_FTA=_b[fta]+_b[l1fta]+_b[l2fta]
				matrix CVM=e(V)
				matrix list CVM
				global cov_fta_l1fta=CVM[2,1]
				global cov_fta_l2fta=CVM[3,1]
				global cov_l1fta_l2fta=CVM[3,2]
				
				global tot_FTA_var=_se[fta]^2+_se[l1fta]^2+_se[l2fta]^2+2*$cov_fta_l1fta+2*$cov_fta_l2fta+2*$cov_l1fta_l2fta
				global tot_FTA_SE=sqrt($tot_FTA_var)
				
				disp $tot_FTA
				disp $tot_FTA_SE
				
			lincom fta+l1fta+l2fta
			
			global tot_FTA_SE=r(se)
				disp $tot_FTA_SE
			
			* export results using outreg 
			if "`fe'"=="none" & `level'==0 {
				outreg2 using "BB_lags_`method'.txt", excel keep($reglist) cttop(`level'd with `fe') addstat(Total FTA effect, $tot_FTA, SE total FTA effect, $tot_FTA_SE) replace
			}
			else {
				outreg2 using "BB_lags_`method'.txt", excel keep($reglist) cttop(`level'd with `fe') addstat(Total FTA effect, $tot_FTA, SE total FTA effect, $tot_FTA_SE) append
			}
						
			* export coefficient estimates
			keep if _n==1
			foreach var in $reglist {
				gen beta_`var'=_b[`var']
			}
			gen hs="`level'-digit"
			gen fe="`fe'"
			gen method="`method'"
			keep beta* hs fe method
			order method hs fe beta*
			save "coeff_`level'_`fe'_`method'", replace
		
		}		
   }
}

log close

****************************************************************
*** DIFFERENT AGGREGATION LEVELS, HETEROGENEOUS COEFFICIENTS ***
************ Generate results to replicate Figure 1 ************
****************************************************************

cd "$resultspath"
cls
capture log close
local logname="BBComtrade_hetcoeff_"+"$S_DATE"
local logname=subinstr("`logname'"," ","",.)
log using "`logname'", replace

*global fetype="exp_year imp_year"
global fetype="exp_year_product imp_year_product imp_exp_product"

	foreach method in ols ppml {
	foreach level in 2 4 {

		foreach regs in "fta l1fta l2fta" {
		
		disp("`method'")
		disp("`level'")
		disp("$fetype")
				
			* create dataset at appropriate aggregation level
			use "$dtapath\BB_allyears_all_zeros_2000obsplus.dta", replace
			cap gen sitc_`level'_digit=substr(commoditycode,4,`level')
			qui drop if year>2000
			if `level'<4 {
				collapse (sum) imports (mean) comlang_off comcol curcol col45 distw EIA WTO_j WTO_i WTO PTA scott_fta fta, by(importeriso3 exporteriso3 sitc_`level'_digit year)
			}
			
			* generate lags and log imports
			if "`method'"=="ols" {
					*keep if imports~=0   // o/w we get more lags than in original B&B (2007) paper!
			}
			bys importeriso3 exporteriso3 sitc_`level'_digit (year):  gen l1fta=fta[_n-1] if year==year[_n-1]+5 | (year==1965 & year==year[_n-1]+3)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year):  gen l2fta=fta[_n-2] if year==year[_n-2]+10 | (year==1970 & year==year[_n-2]+8)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year):  gen f1fta=fta[_n+1] if year==year[_n+1]-5 | (year==1962 & year==year[_n+1]-3)

			qui gen limports=log(imports)
			qui egen exp_year_product=group(exporteriso3 sitc_`level'_digit year )
			qui egen imp_year_product=group(importeriso3 sitc_`level'_digit year)
			qui egen imp_exp_product=group(importeriso3 exporteriso3 sitc_`level'_digit)
			codebook exp_year imp_year 
			qui egen pairid_str=concat(exporteriso3 importeriso3)
			qui encode pairid_str, gen(pairid)
			
			foreach var in fta l1fta l2fta {
				qui gen coeff_`var'=.
				qui gen se_`var'=.
			}

			* now poisson with ppmlhdfe for each SITC code separately
			qui egen sitccount=group(sitc_`level'_digit)
			cap drop temp
			qui egen temp=max(sitccount)
			local countsitc=temp
			disp `countsitc'
			qui drop temp
			save tempfile.dta, replace
					
			forvalues sitccode=1(1)`countsitc' {
				disp `sitccode'
				use tempfile.dta, replace
				qui keep if sitccount==`sitccode'
				qui count if imports~=0
				qui gen noposobs=r(N)
				qui gen noobs=_N
					if "`method'"=="ppml" {
						noisily: capture ppmlhdfe imports `regs' if imports~=., abs($fetype) cluster(pairid) separation(fe mu)
					} 
					if "`method'"=="ols" {
						noisily: capture reghdfe limports `regs' if imports~=., abs($fetype) cluster(pairid)
					} 
					qui gen temp=e(converged)
					qui replace temp=1 if temp==.
					global conv=temp
					drop temp
					*disp $conv
					if $conv==0 {
					exit // (no convergence)
					}
					if $conv==1 {
						foreach var in `regs' {
							capture replace coeff_`var'=_b[`var'] if sitccount==`sitccode'
							capture replace se_`var'=_se[`var'] if sitccount==`sitccode'
						}
					}
					qui bys sitc_`level'_digit: keep if _n==1
					keep sitc_`level'_digit coeff_* se_* noobs noposobs
					qui gen sitcagglevel="`level'-digit"
					qui gen converged=$conv
					order sitcagglevel sitc_`level'_digit
					qui save temp`sitccode'.dta, replace
				}
			}
			
			use temp1.dta, replace
			erase temp1.dta
			forvalues x=2(1)`countsitc' {
				append using temp`x'.dta
				erase temp`x'.dta
			}
			local fetype=subinstr("$fetype"," ","",.)
			save coeffhet_`method'_`fetype'_`level'_digit.dta, replace
		}
	}

erase tempfile.dta
log close

*** ANALYSE COEFFICIENT HETEROGENEITY ***

cls
cd "$resultspath"

* compute sector-level weights for weighting below
foreach level in 2 4 {
disp("`level'")
use "$dtapath\BB_allyears_all_zeros_2000obsplus.dta", replace
cap drop sitc_`level'_digit
gen sitc_`level'_digit=substr(commoditycode,4,`level')
egen totsectortrade=sum(imports), by(sitc_`level'_digit)
egen tottrade=sum(imports)
gen secweight=totsectortrade/tottrade
bys sitc_`level'_digit: keep if _n==1
keep sitc_`level'_digit secweight
save "secweights_`level'_digit.dta", replace
}

* Kernel density of coefficient estimates at 2-digit and 4-digit level

cls
local lb=1e6
local ub=-1e6
foreach method in ols ppml {
foreach level in 2 4 {
use coeffhet_`method'_exp_year_productimp_year_productimp_exp_product_`level'_digit.dta, replace
qui gen xaxis=coeff_fta+coeff_l1fta+coeff_l2fta
su xaxis, meanonly
local lb=min(r(min),`lb')
local ub=max(r(max),`ub')
}
}
clear
set obs 500
local lb=`lb'*1.05
local ub=`ub'*1.05
local delta=(`ub'-`lb')/499
egen xaxis=fill(`lb'(`delta')`ub')
mkmat xaxis, matrix(xaxis)

foreach method in ols ppml {
foreach level in 2 4 {
disp ("`method'")
disp("`level'")
	
* homogeneous coefficient estimates		
use coeff_0_imp_exp_bilateral_`method'.dta, replace // This should be zero!
		global homcoeff=beta_fta+beta_l1fta+beta_l2fta
		disp $homcoeff
		use coeffhet_`method'_exp_year_productimp_year_productimp_exp_product_`level'_digit.dta, replace
bys sitc_`level'_digit: keep if _n==1	
		gen coeff_fta_total=coeff_fta+coeff_l1fta+coeff_l2fta
		label var coeff_fta_total "Total FTA effect (log points), 2-digit"
		sum coeff_fta_total, det
		
		*     
		* trade-weighted average of coefficient estimates
		disp("`level'")
		merge 1:1 sitc_`level'_digit using "secweights_`level'_digit.dta
		drop _merge
		gen temp=coeff_fta_total*secweight
		egen coeff_fta_total_weighted=sum(temp)
		global coeff_fta_total_weighted=coeff_fta_total_weighted
		disp $coeff_fta_total_weighted
		svmat xaxis, names(xaxis)
		su xaxis coeff_fta_total
		di $coeff_fta_total_weighted "   " $homcoeff
		
		if "`method'"=="ppml" kdensity coeff_fta_total, at(xaxis) bwidth(0.2) xline($homcoeff, lwidth(medium)) xline($coeff_fta_total_weighted, lp(dash) lwidth(medium)) title(" ")xtitle("PPML estimation, `level'-digit coefficients", margin(medium)) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) xlabel(-2(1)3) ylabel(0(0.25)1.10, grid) saving(hist_`level'_digit_`method', replace) 
		else kdensity coeff_fta_total, at(xaxis) bwidth(0.2) xline($homcoeff, lwidth(medium)) xline($coeff_fta_total_weighted, lp(dash) lwidth(medium)) title(" ")xtitle("OLS estimation, `level'-digit coefficients", margin(medium)) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white)) xlabel(-2(1)3) ylabel(0(0.25)1.10, grid) saving(hist_`level'_digit_`method', replace) 
	}
}
gr combine hist_2_digit_ols.gph hist_4_digit_ols.gph hist_2_digit_ppml.gph hist_4_digit_ppml.gph  , col(2) row(2) graphregion(fcolor(white) lcolor(white) ifcolor(white) ilcolor(white))

capture for X in any 2 4: for Y in any ols ppml: erase hist_X_digit_Y.gph
graph save "Graph" "Figure 1.gph", replace


************************************************************
*** COMPARE FORECASTS AT DIFFERENT LEVELS OF AGGREGATION ***
************ Replicate results in the Appendix *************
************************************************************
cls
cd "$resultspath"

*global fetype="imp_exp"  // merges in file base on FE: exp_year imp_year
global fetype="all"  // merges in  file based on FE: exp-year-product, imp-year-product, and imp-exp-product

foreach level in 0 2 4  {
	foreach method in ols ppml {
		foreach type in het hom {

		disp("`level'")
		disp("`method'")
		disp("`type'")
		disp("$fetype")

			* create dataset at appropriate aggregation level
use "$dtapath\BB_allyears_all_zeros_2000obsplus.dta", replace
			
			  cap drop  sitc_`level'_digit
			gen sitc_`level'_digit=substr(commoditycode,4,`level')
			replace sitc_`level'_digit="0" if sitc_`level'_digit==""
			drop if year>2000
			if `level'<4 {
				collapse (sum) imports (mean) fta, by(importeriso3 exporteriso3 sitc_`level'_digit year)
			}
					
			* generate lags and log imports
			if "`method'"=="ols" {
				*keep if imports~=0   // o/w we get more lags than in original B&B (2007) paper!
			}
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen l1fta=fta[_n-1] if year==year[_n-1]+5 | (year==1965 & year==year[_n-1]+3)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen l2fta=fta[_n-2] if year==year[_n-2]+10 | (year==1970 & year==year[_n-2]+8)
			bys importeriso3 exporteriso3 sitc_`level'_digit (year): gen f1fta=fta[_n+1] if year==year[_n+1]-5 | (year==1962 & year==year[_n+1]-3)
			gen limports=log(imports)
			
				* merge in coefficient estimates
				
				if "`type'"=="hom" {
				gen method="`method'"
				local fetype="imp_exp_bilateral" // note: if we impose coefficient homogeneity, we also remove the product-dimension from the FE!
				merge n:1 method using "$resultspath\coeff_`level'_`fetype'_`method'.dta"
				drop _merge
							
				for any fta l1fta l2fta: ren beta_X coeff_X
			}
		
			if "`type'"=="het" {
				if `level'~=0 {
					local fetype="$fetype"
					if "`fetype'"=="imp_exp_prod" {                            
						merge n:1 sitc_`level'_digit using "$resultspath\coeffhet_`method'_exp_yearimp_year_`level'_digit.dta"
						
					}
					if "`fetype'"=="all" {
					disp("entering merge loop")
						merge n:1 sitc_`level'_digit using "$resultspath\coeffhet_`method'_exp_year_productimp_year_productimp_exp_product_`level'_digit.dta"
									
				}
					 drop _merge
				
				}
				if `level'==0 {
			disp("entering zero loop")
				gen method="`method'"
				local fetype="$fetype"
				merge n:1 method using "$resultspath\coeff_`level'_`fetype'_`method'.dta"
				drop _merge
				
				for any fta l1fta l2fta: ren beta_X coeff_X

				}
			}
			
			* predict change in tradeflows due to FTAs after 10 years
			cls
			foreach teff in att atu ate {
			*foreach teff in atu {
				disp("`teff'")
				preserve
				if "`teff'"=="att" {
				* if we compute ATT, we must only keep countrypairs with FTAs that are fully phased in (so the FTA dumy and its lags are all 1)
					*bys importeriso3 exporteriso3 sitc_`level'_digit (year): keep if l2fta==1 & l2fta[_n-1]==0  // 
					keep if l2fta==1 & l1fta==1 & fta==1 
				}
				if "`teff'"=="atu" {
				
				* if we compute ATU, we must only keep countrypairs without FTAs (so the FTA dummy and all its lags are =0)
					keep if l2fta==0 & l1fta==0 & fta==0 
				}
				sum fta l1fta l2fta
				*foreach var in fta l1fta l2fta {
				*	replace coeff_`var'=0 if coeff_`var'/se_`var'<1.65  // only use coefficient at 10% significance or more
				*}
				*gen double dlimports=coeff_fta*fta+coeff_l1fta*l1fta+coeff_l2fta*l2fta
				*gen double imports_cf=exp(limports-dlimports)    // still need to add a term reflecting the exponent of the error term?
				*gen double imports_cf=(imports/dimports) if imports~=0
				*sum dlimports, det
				gen double dimports=exp(coeff_fta+coeff_l1fta+coeff_l2fta) if imports~=0
				sum dimports, det
				
				if "`teff'"=="att" {
					gen double imports_cf=(imports/dimports) if imports~=0
					egen double aggimp_cf=sum(imports_cf)
					egen double aggimp=sum(imports)
					gen double te=(aggimp/aggimp_cf)-1
				}
				if "`teff'"=="atu" {
					gen double imports_cf=(imports*dimports) if imports~=0
					egen double aggimp_cf=sum(imports_cf)
					egen double aggimp=sum(imports)
					gen double te=(aggimp_cf/aggimp)-1
				}
								if "`teff'"=="ate" {
				
				* ate is increase over all country pairs, starting from situation without FTA
				* so need to compute: 1) sum of trade over all country pairs without FTAs 
				* and 2) sum of trade over all country pairs with FTAs
				gen double imports_fta=imports if (l2fta==1 & l1fta==1 & fta==1 ) // if countries already have an FTA, use this for the first sum
				replace imports_fta=(imports*dimports) if imports~=0 &(l2fta==0 & l1fta==0 & fta==0) // if countries don't have an FTA, need to compute a counterfactual for the first sum
								
				gen double imports_nofta=imports if  (l2fta==0 & l1fta==0 & fta==0) // for countries without an FTA, simply use current bilateral trade for the second sum
				replace imports_nofta=(imports/dimports) if imports~=0 & (l2fta==1 & l1fta==1 & fta==1) // if countries already have an FTA, need to compute a counterfactual for the second sum
				egen tottrade_fta=sum(imports_fta)
				egen tottrade_nofta=sum(imports_nofta)
				gen double te=(tottrade_fta/tottrade_nofta)-1
				}
				sum 
				
				keep te
				disp "`teff'"
				sum
				
				keep if _n==1
				gen level=`level'
				gen method="`method'"
				gen type="`type'"
				gen tetype="`teff'"
				local fetype="$fetype"
				save `teff'_`method'_`level'_`type'_`fetype'.dta, replace
				restore
				

			}
		}
	}
}

* produce treatment effects table
global fetype="all"
cls

foreach level in 0 2 4  {
	foreach method in ols ppml {
		foreach type in het hom {
			foreach teff in att atu ate {
				local fetype="$fetype"
				if `level'==0 & "`method'"=="ols" & "`type'"=="het" & "`teff'"=="att" use `teff'_`method'_`level'_`type'_`fetype'.dta, replace
				else append using `teff'_`method'_`level'_`type'_`fetype'.dta			
			                           }
		                         }
	                            }
                         }
cd "$resultspath"						 
save "Treatment Effects.dta", replace

	
	
	





